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ABSTRACT 


A method has been developed for real-time stability margin measurement calculations. 
The method relies on a tailored-forced excitation targeted to a specific frequency range. 
Computation of the frequency response is matched to the specific frequencies contained in the 
excitation. A recursive Fourier transformation is used to make the method compatible with real-time 
calculation. The method was incorporated into the X-38 nonlinear simulation and applied to an 
X-38 robustness test. X-38 stability margins were calculated for different variations in aerodynamic 
and mass properties over the vehicle flight trajectory. The new method showed results comparable 
to more traditional stability analysis techniques, and at the same time, this new method provided 
coverage that is more complete and increased efficiency. 
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total time (window length) to be processed rounded to nearest discrete time step, sec 

time series of experimental excitation input applied to control effector actuator, deg 

k lh discrete frequency in excitation, rad/sec 

desired minimum frequency of excitation, cycles/sec 

series of data in time domain 

complex number of data transformed to frequency domain 
time separation between points in measure data, sec 
circumference of circle divided by two times its radius 
Schroeder phase shift to be applied at k th frequency, rad 
frequency, rad/sec 

frequency of high-pass filter, rad/sec 

Functions 

take imaginary part of complex number 
base ten logarithm of number 
take real part of complex number 
arctangent of angle 


INTRODUCTION 


The conventional method for testing robustness of a space vehicle is to apply a Monte Carlo 
analysis. This involves varying some parameters, such as the aerodynamic, environmental, sensor, 
and/or mass properties of the vehicle; flying the vehicle trajectory; and assessing the percentage of 
success or failure. The success or failure criterion is derived from a general assessment that the 
mission goal was achieved. For successful instances, however, there is little indication of how close 
to failure each particular instance was. Adding a stability margin calculation to Monte Carlo 
analysis would provide a better indication of the true robustness of the vehicle. 

The X-38 project developed a robustness test that combined some features of Monte Carlo 
analysis with conventional stability analysis. To accomplish this analysis combination, a set of test 
instances was selected with specific variations in aerodynamic and mass properties. Stability margin 
requirements were placed on each uncertainty instance. The stability margins were calculated using 
either a high-fidelity linear model or by frequency sweeps with the nonlinear simulation in a 
five-degree-of-freedom mode. This analysis was relatively time-consuming and relied on validated 
linear models or long duration frequency sweeps. In addition, since only two points along the 
trajectory were analyzed, there was no guarantee that the worst-case flight condition was addressed. 
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Therefore, to increase efficiency and cover more of the X-38 flight conditions, a new approach 
was developed that incorporates a real-time stability (RTS) margin calculation into the full 
six-degree-of-freedom X-38 nonlinear simulation. Results using the new method were compared to 
the original X-38 robustness test results. This report discusses this comparison, the improvements 
to the real-time stability margin measurement algorithm, and the benefits and weaknesses of this 
technique. 


BACKGROUND 


Linear stability analysis has been used extensively to test for vehicle robustness. The X-29 
program used fast Fourier transformation (FFT) techniques to measure stability margins in flight 
(ref. 1). The pilot-generated frequency sweeps used for this in-flight measurement required long 
periods of stabilized flight (approximately one min). 

Space vehicles such as the X-38 or X-34 are typically in a state of rapidly accelerating or 
decelerating flight. For vehicles with rapidly changing flight conditions, a tailored-forced excitation 
method was developed (ref. 2) to reduce the requirement for long periods of stabilized flight. 

It was desired to incorporate the stability margin measurement into the Monte Carlo testing that 
is commonly used to test the robustness of space vehicles. By including a stability margin 
measurement into the Monte Carlo analysis, instances that may have survived but with very low 
margin, could be identified. In addition, sensitivity of the stability margin with respect to different 
parameters could be assessed. Conventional FFT methods combined with the large number of 
instances in a typical Monte Carlo run results in an impractical computational burden. Fortunately, 
a recursive Fourier transformation algorithm (ref. 3) can be used to dramatically reduce the 
computational requirements of the stability margin calculation. The subject of this report is this 
real-time (online) algorithm and its application to the X-38 Monte Carlo analysis. 

X-38 PROGRAM BACKGROUND 


The goal of the X-38 program was to demonstrate the technology required for the International 
Space Station’s emergency crew return vehicle (CRV) (ref. 4). Based on the X-24A lifting body, 
five prototype vehicles were developed to demonstrate the CRV concept. These unmanned vehicles 
operate autonomously and use a steerable parachute, called a parafoil, for final descent and landing. 
Of the five X-38 vehicles, four were to be used only for atmospheric flight-testing and one was to 
be tested from space through reentry to landing. As these vehicles had no propulsion systems, they 
were to be carried to altitude under the wing of a B-52 airplane and released. 

Vehicle 131R is one of the atmospheric flight test vehicles (fig. 1). The robustness analysis 
presented in this report was performed using the nonlinear simulation of vehicle 131R during free 
flight. Free flight (ref. 4) is the portion of flight after release from the B-52 airplane but before 
deployment of the parafoil. During free flight, the X-38 vehicle uses a dynamic inversion flight 
control system (ref. 5). A nominal free-flight trajectory lasts for 55 sec starting at Mach 0.75 at an 
altitude of 45,000 ft and ends at Mach 0.74 at an altitude of 28,000 ft (fig. 2). 
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As part of preflight certification, the X-38 vehicle 13 1R was required to pass a robustness test. 
A set of 20 test instances (20 uncertainty instances) was selected with specific parametric variations 
in aerodynamic and mass properties (table 1). Based on previous experience with lifting bodies and 
the space shuttle, values for three-sigma variation of these parameters were defined. The values in 
table 1 show the percent of this three-sigma change applied to each parameter. For each uncertainty 
instance, the vehicle stability margin was tested at high and low dynamic pressure flight conditions. 
The low dynamic pressure condition was defined as the point 8 sec into the free-flight trajectory. 
The high dynamic pressure condition was defined as 47 sec into the free-flight trajectory. The 
control system passed the robustness test if the margins at these two conditions were greater than 
4-dB gain margin and 20-deg phase margin for all of the uncertainty test instances. The nominal 
instance required 6-dB gain margin and 45 -deg phase margin. 

Johnson Space Center performed this robustness analysis with a high-fidelity linear model. 
Dryden Flight Research Center (DFRC) independently verified the results. The DFRC analysis 
used an FFT approach with a nonlinear simulation. The simulation was limited to 
five-degree-of-freedom to maintain the altitude and Mach number during the relatively long 
duration frequency sweep excitation. Results from this FFT-based approach provide the basis for 
comparison of the new real-time technique. 



Figure 1. X-38 vehicle 131R. 
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Figure 2. X-38 nominal free-flight trajectory. 
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Table 1. X-38 vehicle 13 1R uncertainty test instances. 


Case name 

la 

lb 

If 

2a 

2b 

2f 

3a 

3b 

3d 

3e 

Yaw due to beta 

-66.7 

-66.7 

-66.7 

-66.7 

-66.7 

-66.7 

+66.7 

+66.7 

+66.7 

+66.7 

Roll due to beta 

-66.7 

-66.7 

+66.7 

+66.7 

+66.7 

+66.7 

+66.7 

+66.7 

+66.7 

+66.7 

Roll due to aileron 

-66.7 

-66.7 

-66.7 

-66.7 

-66.7 

-66.7 

-66.7 

+66.7 

-66.7 

-66.7 

Yaw due to aileron 

-66.7 

-66.7 

-66.7 

-66.7 

-66.7 

-66.7 

-66.7 

+66.7 

-66.7 

-66.7 

Roll due to rudder 

-66.7 

+66.7 

+66.7 

+66.7 

-66.7 

-66.7 

+66.7 

+66.7 

+66.7 

-33.3 

Yaw due to rudder 

+66.7 

-66.7 

-66.7 

-33.3 

+66.7 

+66.7 

-33.3 

-33.3 

-33.3 

+66.7 

Side force due to beta 

-66.7 

-66.7 

-66.7 

-66.7 

-66.7 

-66.7 

+66.7 

+66.7 

+66.7 

+66.7 

Side force due to rudder 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

Side force due to aileron 

-66.7 

-66.7 

-66.7 

-66.7 

-66.7 

-66.7 

-66.7 

+66.7 

-66.7 

-66.7 

Pitching moment 

-66.7 

+66.7 

+66.7 

-66.7 

+66.7 

+66.7 

-66.7 

+66.7 

+66.7 

-66.7 

Pitching moment due to elevon 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

Rolling moment 

+66.7 

+66.7 

+66.7 

+66.7 

+66.7 

+66.7 

+66.7 

+66.7 

+66.7 

+66.7 

Yawing moment 

+66.7 

+66.7 

+66.7 

+66.7 

+66.7 

+66.7 

+66.7 

+66.7 

+66.7 

+66.7 

Longitudinal CG change 

forward 

aft 

aft 

forward 

aft 

aft 

forward 

aft 

aft 

forward 

Lateral CG change 

right 

right 

right 

right 

right 

right 

right 

right 

right 

right 

Vertical CG change 

up 

down 

up 

up 

down 

up 

up 

down 

down 

up 


Table 1. Concluded. 


Case name 

4a 

4b 

4d 

5a 

5b 

6a 

6b 

7a 

8a 

8b 

Yaw due to beta 

+66.7 

+66.7 

+66.7 

-66.7 

-66.7 

-66.7 

-66.7 

-66.7 

-81.0 

-66.7 

Roll due to beta 

-66.7 

-66.7 

-66.7 

-66.7 

-66.7 

+66.7 

+66.7 

-66.7 

0.0 

+66.7 

Roll due to aileron 

+66.7 

+66.7 

+66.7 

+66.7 

+66.7 

-66.7 

-66.7 

+66.7 

0.0 

0.0 

Yaw due to aileron 

+66.7 

+66.7 

+66.7 

+66.7 

-66.7 

+66.7 

+66.7 

-66.7 

0.0 

0.0 

Roll due to rudder 

+66.7 

-66.7 

-33.3 

+66.7 

+66.7 

-66.7 

-66.7 

-66.7 

0.0 

+66.7 

Yaw due to rudder 

-33.3 

+33.3 

+66.7 

-66.7 

-66.7 

+66.7 

+66.7 

+66.7 

0.0 

-66.7 

Side force due to beta 

+66.7 

+66.7 

-66.7 

-66.7 

-66.7 

-66.7 

-66.7 

-66.7 

+81.0 

+66.7 

Side force due to rudder 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

Side force due to aileron 

+66.7 

+66.7 

+66.7 

+66.7 

-66.7 

+66.7 

+66.7 

-66.7 

0.0 

0.0 

Pitching moment 

-66.7 

+66.7 

+66.7 

-66.7 

+66.7 

-66.7 

+66.7 

-66.7 

0.0 

0.0 

Pitching moment due to elevon 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

Rolling moment 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

Yawing moment 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

Longitudinal CG change 

forward 

aft 

forward 

forward 

aft 

forward 

aft 

aft 

none 

none 

Lateral CG change 

right 

right 

right 

right 

right 

right 

right 

right 

none 

none 

Vertical CG change 

up 

down 

down 

up 

down 

up 

down 

up 

none 

none 
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REAL-TIME STABILITY MARGIN CALCULATION APPROACH 


Reference 2 outlines a test technique to allow for measuring the stability of a vehicle or flight 
control system while minimizing the duration of the required input excitation. This report addresses 
how that technique was modified to be more compatible with real-time (online) processing. 
The major change was to incorporate a recursive Fourier transform algorithm to significantly reduce 
the computational requirements. 

The real-time stability method uses a forced excitation applied to a control effector command. 
Measured time segments (windows) of the appropriate input and output variables are transformed 
to the frequency domain, and the open-loop frequency response is formulated. The resulting 
frequency response is used to calculate stability margins of the vehicle or flight control system. A 
time history of the stability margin can be obtained by repeating this calculation at each new time 
step. The time segment is moved by one time step for each new calculation, providing a sliding 
window of data. 


Excitation Signal 

The excitation signal used is a summation of cosine waves of specially selected frequencies. 
The frequency range selected is application specific and is further defined by the expected 
frequency at which the lowest stability margin occurs (crossover frequency). By carefully matching 
specific frequencies and the length of time segment (window length) being analyzed, problems 
associated with frequency leakage (ref. 2) are avoided, and good quality stability margin 
measurements are obtained. 


Window Length 


The lowest desired frequency and the number of cycles of that frequency determine the length 
of the time segment (window length) used for analysis. The stability margin calculation is a result 
of the dynamic response over the window length. A longer window length adds an effective delay of 
half the window length to the “real-time” calculation. Shorter window lengths are also desired for 
minimizing flight condition change during margin estimation. Window length may need to be 
increased, however, for better frequency resolution or to improve results in the presence of noisy 
data. In these cases, using multiple cycles at the lowest frequency lengthens process time. 


The process time is computed from the number of cycles at the desired minimum frequency 
using: 


T 


win 


1 

W m i n ( c Y cles/sec ) 


X N cycles 


( 1 ) 


This process time is rounded to the nearest discrete time step (At), which is usually the update 
rate of the flight control laws. The starting frequency is then adjusted to fit the actual window 
length. Higher frequencies whose cycle lengths fit evenly in the window defined above are selected 
using: 


W^rad/sec) = 2tzx 



win 


, where k is any positive integer. 


( 2 ) 
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Phase Shifting 


The top plot in figure 3 shows a typical excitation generated by a sum of cosines with no phase 
shifting of the individual frequencies. As suggested by Schroeder (ref. 6), the phase of each of the 
cosine waves was shifted before summing to reduce the peak factor of the excitation signal while 
maintaining the frequency content. The amount of phase shifting to apply to each of the discrete 
frequencies is computed from: 





The total excitation input is then obtained using the following equation: 


(3) 


U(t) = [£cos(Wk^t + ^ k )VN f (4) 

k = 1 

The division by the total number of frequency components in excitation input (Nj ) is used to 
normalize (scale) the amplitude of the resultant signal. The bottom plot on figure 3 shows an 
example of an excitation signal with Schroeder-phase shifting. 




Figure 3. Effect of Schroeder-phase shifting on excitation signal. 
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Data Filtering 


An important step in calculating a transfer function using FFT techniques is to de-trend the 
input and output time histories. The simplest way to perform this de-trend is to subtract the average 
of all the data points from each data point. This step can significantly improve the results when the 
data has a nonzero steady-state value (i.e., an elevator excitation with a nominal trim value of 5 deg). 
This method works well for a single frequency response calculation over a single time period. The 
objective is to calculate the frequency response at each time step using a moving window of data. To 
calculate an average at each time step would require storing many data points and computing an 
average at each time step. 

A more efficient approach is to apply a high-pass filter to both the input and output data stream 
before doing the transformation to the frequency domain. The high-pass filter effectively removes 
the steady-state component of the time history data. Since the same filter is applied to both input 
and output, the transfer function relationship (ratio of output due to input) is preserved. The 
following simple high-pass filter was chosen: 


to HP S 
S “T CO j-^p 


( 5 ) 


A value of 0 ) Hp = 6 rad/sec seemed to work well. For u> Hp = 6 rad/sec, frequencies above 
1.0 rad/sec are amplified and frequencies below 1.0 rad/sec are attenuated. This value worked for 
the X-38 vehicle where the lowest crossover frequencies are around 3 rad/sec, but this value might 
need adjustment for other vehicles with very low crossover frequencies. 


Recursive Fourier Transformation 

If the stability margin is calculated at each time step, using an FFT algorithm quickly becomes 
computationally prohibitive. Even if a Chirp-Z algorithm (ref. 2) is used with a limited number of 
frequencies, the calculations generally exceed the frame time. One solution is to run the calculation 
in the background, computing a margin at a slower update rate using leftover time spread over a 
number of computer timeframes. 

Alternatively, a more efficient algorithm, such as the recursive Fourier transform, can be used 
for transformation to the frequency domain (ref. 3). This algorithm was used when performing the 
analysis on the X-38 vehicle. The increased algorithm efficiency allowed for stability margin 
calculation at each time step of the X-38 simulation. 

The definition of the finite Fourier transform is: 

V 

X(co ) = j x(t) e-jnt dt (6) 
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Thus, the transformation from time domain to frequency domain is accomplished by integrating 
over the time segment of interest. When the sample time ( At ) is short compared to the period of the 
frequencies of interest, a simple Euler integration approximation can be used: 


m 

X(j(0) = Afx L (7) 

i = 1 

Integration is required for each frequency of interest. With this transform, there is no restriction 
on which frequencies can be used, so the excitation frequencies can be used exactly. 

To save computation time, it is not necessary to repeat this summation for each time segment. 
Instead, the discrete data sample points are stored in a buffer. At each new time step, the new data 
sample is added and the oldest data point is subtracted. In this manner, the process time is 
significantly reduced. This computational simplification comes at the cost of the memory required 
to save a potentially large number of data points in a buffer. 


Computing Frequency Response 

After the input and output have been transformed to the frequency domain, the two complex 
vectors are divided to create the open-loop transfer function: 


TF(j co) 


X , , 

output 

V 

input 


(8) 


The elements of the complex vectors are then converted from complex number form to 
magnitude and phase form using the following: 


phase = tan' 1 WFC/m)) 
Re(TF(j'co)) 


(9) 


magclB = 20 x log 10 ^ Re (rF(jco)) 2 + lm(TF(ja>)) 2 ) ( 10 ) 

The magnitude and phase of the transfer function can then be plotted against the frequency 
vector to form points on a Bode plot. Bode plot points calculated by the real-time stability method 
have compared well against Bode plots from the frequency sweep method as shown in figure 4. 
This result is a typical match for a linear system with good excitation. 
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Figure 4. Comparison of frequency responses obtained from real-time stability and frequency 
sweep methods. 


Stability Margin Calculation 

Stability margins are calculated from the open-loop frequency response data using simple linear 
interpolation. Each point in the magnitude vector is compared with neighboring points to find two 
points that enclose the 0-dB crossover point. Linear interpolation is then used to estimate the gain 
crossover frequency and phase margin. To find the gain margin, points surrounding the -180-deg 
phase crossover are located. Linear interpolation is again used to estimate the actual phase crossover 
frequency and gain margin. 

Errors introduced by linear interpolation are directly related to the frequency resolution of the 
frequency response points. The resolution can be improved by increasing the length of the process 
window. To achieve better resolution, either choose a lower value for the lowest desired frequency 
or increase the number of cycles at the lowest frequency. 
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APPLICATION TO X-38 ELEVATOR LOOP 


The gain crossover for pitch axis of the X-38 vehicle typically occurs at around 3.5 rad/sec or 
higher. To obtain good margin calculations and sufficient frequency resolution, the following 
values were chosen: 


W ■ = 0.477 cycles/sec (3.0 rad/sec) 


(ID 


^ cycles ^ 


( 12 ) 


The X-38 control system calculation rate is 25 times per sec. These parameters resulted in a time 
segment window length of 4.20 sec (105 points). An excitation signal was generated using the 
first nine frequencies with wavelengths that matched the time segment window. Table 2 shows these 
frequencies. 


Table 2. Frequency content 
of X-38 excitation. 


k 

W k , rad/sec 

1 

2.9920 

2 

4.4880 

3 

5.9840 

4 

7.4800 

5 

8.9760 

6 

10.4720 

7 

11.9680 

8 

13.4640 

9 

14.9600 


Figure 2 shows a time history of a nominal X-38 trajectory. The vehicle is released from under 
the wing of a B-52 airplane at Mach 0.75 and an altitude of 45,000 ft. After an initial transient 
resulting from aerodynamic interaction with the B-52 airplane, the vehicle stabilizes at a commanded 
16-deg angle of attack. At around 25 sec, the commanded angle of attack is changed to 18 deg. 
The vehicle flies for 55 sec of free flight before the parafoil deployment sequence is initiated. 
The X-38 stability robustness test was performed over this 55-sec period of free flight. 
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The dashed line on figure 2 also shows the same free-flight trajectory with the previously 
described excitation applied symmetrically to both elevators. The amplitude of the excitation was 
scaled to +0.5 deg magnitude. The applied excitation provided small perturbations to pitch rate and 
angle of attack but did not significantly affect the vehicle trajectory. 

The X-38 vehicle was flown on this trajectory with 20 different aerodynamic uncertainty 
variations (table 1). Figure 5 shows the resulting elevator loop phase and gain margins with their 
associated crossover frequencies. After approximately 10 sec, the calculated margins settle out to 
consistent reasonable results. Transients associated with release from the B-52 airplane corrupt 
results before 10 sec. 

It should be noted that the stability margin is a calculation based on a time segment of data. In 
this case, the time segment (window) length is 4.2 sec. The time histories in figure 5 show a margin 
that is the result of the previous 4.2 sec of data. The margin is effectively delayed by half the 
window length. This provides another reason to use a shorter window length. 



Figure 5. X-38 real-time stability margin calculation for elevator loop. 
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This figure 5 data show how stability changes with flight condition (and corresponding control 
system gain changes) and with varying uncertainties. The worst-case flight condition and 
uncertainty set can be quickly determined. For example, the uncertainty sets show two groupings 
with respect to gain margin. A look at the uncertainties shows that the group with the lower gain 
margins is associated with instances where pitching moment was increased by 66.7 percent. 

To provide validation of these results, a comparison was made between the real-time 
stability method and more conventional frequency response methods using FFT. The flight 
condition obtained at 47 sec (high dynamic pressure) along the trajectory was chosen for this 
comparison. At 47 sec, the flight condition was “frozen” by putting the simulation into a 
five-degree-of-freedom mode. In this mode, the altitude and velocity states are not integrated. 
A relatively long duration (60 sec) frequency sweep was applied to the symmetric elevator 
command. Fourier transformation of the resulting response provided a calculation of the open-loop 
frequency response and the resulting stability margins. Figure 6 shows a comparison of the 
margins found by the real-time stability and the frequency sweep methods. In general, gain margins 
compared to within 0.15 dB and the phase margin was within 3 deg. Some of the phase margin 
difference is a result of interpolation errors caused by the lower frequency resolution of the 
RTS data. 
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Figure 6. Elevator loop stability margins. 
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The simulation run time required to obtain the RTS results was much less than the time required 
to obtain the frequency sweep data. The RTS results can be obtained by running the simulation 
up to the time point of interest. For the frequency response results, the simulation was run up to the 
time point of interest and then run an additional 60 sec for the frequency sweep. In this case, the 
RTS results are for the entire trajectory, and the frequency sweep result is for only one time point. 
To obtain frequency sweep results for each time point would require 60 sec of simulation run time 
for each time point, plus the time to run the simulation until each point. 

APPLICATION TO X-38 RUDDER LOOP 

The stability margins associated with the rudder control loop were calculated using the same 
excitation as used for the elevator. Figure 7 shows the margins and crossovers as a function of time 
for the 20 uncertainty instances. The results for the rudder loop show less consistency than for the 
elevator loop. Some of the instances show significant variation with time. The expectation is that 
even though the flight condition might be changing, the vehicle dynamics should provide a relatively 
continuous and smooth change in behavior. Normally the stability margins should vary in a 
continuous and smooth manner. 
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Figure 7. X-38 real-time stability margin calculation for rudder loop. 
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Closer examination of these results shows that relatively small variations in the frequency 
response caused these larger variations in gain margin. Figure 8 shows the rudder loop frequency 
response calculated (at each time step) over a two-sec period. This plot shows how the calculated 
frequency responses lie within a band. At the crossover frequency (8 rad/sec), the variations 
in phase angle were approximately 5 deg (approximately 10 msec). These fairly small variations in 
phase angle caused variations in the crossover frequency of approximately 1 rad/sec. Although the 
variations in the gain were fairly small (within 1 dB), the steepness of the curve results in gain 
margin variations of approximately 2 dB. 



Figure 8. RTS results over two sec with large gain margin variations. 


This erratic behavior shown in the rudder margin measurements is undesirable. For large Monte 
Carlo runs, regions of consistently low margin may be overlooked because attention is diverted to 
spurious peaks. Filtering the data can provide a clearer overall picture. The data from figure 7 was 
passed through the following filter: 


1 

S+ 1 


( 13 ) 
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Figure 9 shows the filtered margins and crossovers for the rudder loop. The filtered data 
provides a better means of looking for overall trends or for picking global margin minimums. The 
fact that the data were erratic and required filtering, however, may be a clue that there is some 
uncertainty in the margin measurement. For a well-behaved system, the system margins should 
change smoothly and relatively slowly with changes in conditions (i.e., dynamic pressure, angle of 
attack, etc.). Abrupt changes in margin should be investigated as potential problem areas. 






Figure 9. X-38 real-time stability margin calculation for rudder loop with filtering. 


An estimation of the variation of the results can be obtained by simply calculating the difference 
between the filtered and unfiltered results. 

margin error estimation = abs (margin - filtered (margin)) (14) 

crossover error estimation = abs (crossover - filtered (crossover)) (15) 
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The mean and maximum error estimation can be used as part of the statistical analysis of Monte 
Carlo results. The mean of the error estimation can identify instances with less consistent results. 
Figure 10 shows the error estimation values for the rudder loop instances. These data show that 
the phase margin estimation has an error range of approximately ±1.5 deg, and the gain margin is 
within approximately ±1.0 dB. Figure 11 shows that the elevator loop gain margin errors were 
within ±0.25 dB; and the phase margin errors were mostly lower, with the exception of the peak 
errors of approximately 3 deg. These peaks (at roughly 25, 30, and 35 sec) correspond to changes 
in angle of attack and occur as the vehicle was in transition from 16-deg through 18-deg angle 
of attack. 
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Figure 10. X-38 estimated errors for rudder loop. 
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Figure 11. X-38 estimated errors for elevator loop. 
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SIMULTANEOUS MEASUREMENT OF ELEVATOR AND 

RUDDER MARGINS 


To further increase efficiency, the elevator and rudder stability margins can be calculated 
simultaneously. This calculation is achieved by dividing the frequencies between the elevator and 
the rudder. The elevator gain and phase crossovers occur in the range from 3.5 through 4.5 rad/sec 
and 9 through 9.5 rad/sec respectively. The rudder gain and phase crossovers occur in the range 
from 3 through 6 rad/sec and 6.5 through 10 rad/sec respectively. The same discrete frequencies 
used in the previously-mentioned example were divided as shown in Table 3. 


Table 3. Frequency content of excitation used for 


simultaneous inputs. 


k 

W k , rad/sec 

Control effector 

1 

2.9920 

Elevator 

2 

4.4880 

Rudder 

3 

5.9840 

Elevator 

4 

7.4800 

Rudder 

5 

8.9760 

Elevator 

6 

10.4720 

Rudder 

7 

11.9680 

Elevator 

8 

13.4640 

Rudder 

9 

14.9600 

Elevator 


The nominal trajectory (with no uncertainties) was flown with this excitation. Figure 12 shows 
the resulting stability margins obtained with simultaneous and separate excitations for the nominal 
instance. To better illustrate the differences, the data shown has not been filtered as was suggested 
in the previous section. Using simultaneous excitation produced results within 3 deg of phase 
margin and 0.5 dB of gain margin. Most of this difference is caused by interpolation with coarser 
data. Increasing the lowest frequency analyzed or the number of cycles of the lowest frequency 
could achieve better data resolution. 
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Figure 12. X-38 real-time stability margin calculation with simultaneous inputs. 
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INSTANCES WITH REDUCED STABILITY 


The X-38 results were obtained with a mature control system that had good robustness 
characteristics. To show the RTS method with low stability, the rudder loop gain was varied from 
-6 dB through +8 dB. Figure 13 shows that the RTS method demonstrates good results even when 
neutral stability is reached. As the neutral stability point is reached, some variation in the results is 
apparent. At this point, the vehicle is reaching position saturation on the rudders, so some 
degradation in the linear relationship is expected. These data, however, show that the RTS method 
provides a good indication of an approaching instability. 
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Figure 13. X-38 real-time stability margin calculation with variation in rudder loop gain. 


22 


CONCLUDING REMARKS 


A real-time stability method was developed for stability analysis over a vehicle trajectory. The 
method is efficient and compatible with online calculations as well as for use with Monte Carlo 
analysis. The discovery was made that filtering the results and calculating error estimation provided 
measurements that were suitable for use in Monte Carlo analysis. The filtered results allowed a 
better means of obtaining a global minimum margin. The error estimation provided a means for 
flagging instances where the data were inconsistent. 

Results from the real-time stability method compared well with results using frequency sweeps. 
Interpolation errors induced by coarser data caused the largest source of error. Simultaneously 
calculating the elevator and rudder margins can obtain a further increase in efficiency. Results 
obtained using simultaneous inputs compared well with results using separate inputs. 

The real-time stability method works well for instances with reduced stability. The rudder loop 
gain was varied by as much as +8 dB. The real-time stability method produced a good calculated 
stability margin all the way to neutral stability. These data demonstrate that the real-time stability 
method can provide a good indication of an approaching instability. 

The real-time stability method shows the potential to identify deficient instances in large Monte 
Carlo runs. This method is also compatible with real-time (online) applications. 


Dryden Flight Research Center 

National Aeronautics and Space Administration 

Edwards, California, July 15, 2004 
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